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ABSTRACT 

Between 2004 and 2009 a sample of 28 X-ray selected high- and intermediate-frequency peaked 
blazars with a X-ray flux larger than 2 /z Jy at IkeV in the redshift range from 0.018 to 0.361 
was observed with the MAGIC telescope at energies above 100 GeV. Seven among them were 
detected and the results of these observations are discussed elsewhere. Here we concentrate on 
the remaining 21 blazars which were not detected during this observation campaign and present 
the 3sigma (99.7%) confidence upper limits on their flux. The individual flux upper limits lie 
between 1.6 % and 13.6 % of the integral flux from the Crab Nebula. Applying a stacking method 
to the sample of non-detections with a total of 394.1 hours exposure time, we find evidence for an 
excess with a cumulative significance of 4.9 standard deviations. It is not dominated by individual 
objects or flares, but increases linearly with the observation time as for a constant source with 
an integral flux level of ~1.5% of that observed from the Crab Nebula above 150 GeV. 

Subject headings: BL Lacertae objects: general — gamma rays: observations 



1. INTRODUCTION 

MAGIC (Major Atmospheric Gamma-ray 
Imaging Cherenkov) is currently a system of two 
17 m telescopes located atop the Roque de los 
Muchachos on the Canary Island of La Palma at 



a IFAE, Edifici Cn., Campus UAB, E-08193 Bellaterra, 
Spain 

b INAF National Institute for Astrophysics, 1-00136 
Rome, Italy 

c Universita di Siena, and INFN Pisa, 1-53100 Siena, 
Italy 

d Technische Universitat Dortmund, D-44221 Dortmund, 
Germany 

c Universitat Autonoma de Barcelona, E-08193 Bel- 
laterra, Spain 

f Universidad Complutense, E-28040 Madrid, Spain 
STJniversita di Padova and INFN, 1-35131 Padova, Italy 
h Inst. de Astrofi'sica de Canarias, E-38200 La Laguna, 
Tenerife, Spain 

'University of Lodz, PL-90236 Lodz, Poland 
JTuorla Observatory, University of Turku, FI-21500 Pi- 
ikkio, Finland 

k Deutsches Elektronen-Synchrotron (DESY), D-15738 
Zcuthen, Germany 

'ETH Zurich, CH-8093 Switzerland 
m Max-Planck-Institut fiir Physik, D-80805 Miinchen, 
Germany 

"Universitat de Barcelona (ICC/IEEC), E-08028 
Barcelona, Spain 

°Universitat Wiirzburg, D-97074 Wiirzburg, Germany 

p Depto. de Astrofisica, Universidad, E-38206 La La- 
guna, Tenerife, Spain 

CUniversita di Udine, and INFN Trieste, 1-33100 Udine, 
Italy 

Institut de Ciencies de l'Espai (IEEC-CSIC), E-08193 
Bellaterra, Spain 

s Inst. de Astrofi'sica de Andaluci'a (CSIC), E-18080 
Granada, Spain 

'Croatian MAGIC Consortium, Institute R. Boskovic, 
University of Rijeka and University of Split, HR-10000 Za- 
greb, Croatia 

"University of California, Davis, CA-95616-8677, USA 
v Inst. for Nucl. Research and Nucl. Energy, BG-1784 

Sofia, Bulgaria 

w INAF/Osservatorio Astronomico and INFN, 1-34143 

Trieste, Italy 

X ICREA, E-08010 Barcelona, Spain 

CUniversita di Pisa, and INFN Pisa, 1-56126 Pisa, Italy 
supported by INFN Padova 

now at: Centro de Investigaciones Energeticas, 
Medioambientales y Tecnologicas 

***now at: Max-Planck-Institut fiir Kernphysik, D-69029 
Heidelberg, Germany 

+ correspondence: hoehne@astro.uni-wuerzburg.de 



2200 m a. s. 1. The observations refered to in this 
study were obtained during the years 2004 - 2009 
when MAGIC was still a single-dish telescope. Its 
234 m 2 tessellated parabolic mirror allows observa- 
tions of VHE ( Fery High Energy) 7-rays between 
-50GeV and lOTeV. 

One key goal of the MAGIC telescope project 
is to determine the properties of extragalactic 
VHE sources, among which the high-frequency 
peaked BL Lacertae objects are the most numer- 
ous. Blazars are a subclass of radio-loud AGN and 
belong to the most extreme and powerful objects 
in the universe. They are characterized by a non- 
thermal broad-band continuum emission which is 
highly variab le on time scales from years down 



to min utes ([Albert et al.l l2007t lAharonian et al 
2007ah . 



The spectral energy distribution (SED) of 
blazars is characterized by two bumps in a v¥ v 
representation. The first component peaks at 
energies between IR and hard X-rays, and is 
assumed to originate from leptonic synchrotron 
radiation. The maximum of the second peak 
lies in the 7-ray energy regime. The origin of 
this peak can be explained by different and par- 
tially concurring models either rel ying on inverse 



Compton scattering of electrons ([M arasch i et al 



1992; Dermer fc Schlickciser 1993; Sikora et al 



19941 ) or proposing hadro n ic interactions inside 



the j et ( Mannheim 19931 ; iMuecke fc Protheroe 
20011 ). In the case the synchrotron peak oc- 
cu rs at energies above ~ 10 16 5 Hz, (according 
to Nieppola et al. 2006() these blazars are called 
HBLs (high-frequency peaked BL Lacertae ob- 
jects) and for peak energies of ~ 10 14 5-16 5 Hz 
IBLs (intermediate BL Lacertae objects). 

As of April 2010, altogether 29 blazars were es- 
tablished as VHE sources (24 of them HBLs in- 
cluding M87 as 'misaligned' blazar j]], compared to 
six HBLs, when the MAGIC telescope began its 
regular observations in December 2004. The sam- 
ple presented here comprises 21 X-ray selected ob- 
jects which were not detected in the VHE regime 
prior to the MAGIC observations. Nine of the 
objects were already observed between December 
2004 and February 2006 and the up per limits of 
these observations are reported in lAlbert et al 



1 cf . http: / / wwwmagic.mppmu.mpg.de / ~rwagner /sources / 
for an up-to-date list. 
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(|2008al ) . As there have been improvements within 
the MAGIC analysis, the data of these objects 
were re-analyzed and the new results are presented 
in this work. Since no significant detection was at- 
tained, upper limits on a 3 ex (99.7%) confidence 
level will be presented. 

None of the observed sources showed any vari- 
ability on diurnal timescales in the VHE regime. 
Assuming a positive detection in the case of a flar- 
ing state, the observations presented here provide 
a means of investigating the baseline emission of 
these objects. Therefore, a stacking method ap- 
plied to the blazar sample can reveal such an emis- 
sion below the sensitivity limit for each individua l 
object. Together with VERITAS ((Benbowl [2009h 
this is the second stacking analysis which turns 
out to be successful in the VHE 7-ray regime. 
Former experiments like HEGRA failed in de- 
tecting a significant signal in a stacking analysis 
due to their limited s ensitivity (cf. for instance 
Mannheim et al.lll996h . 

In Section [5] the selection criteria for the objects 
will be presented. The observations and the data 
analysis technique are described in Section [3] The 
analysis results are shown in Sectional Finally, a 
discussion of the results and inherent implications 
can be found in Section [5] 

2. BLAZAR SAMPLE 

We selected blaza rs fr om the compilations from 
Donato et al. (2001) and lCostamante fc Ghisellinil 
(|2002n . Additionally, some objects were chosen 
based on the synchr otron peak luminosity from 
Nieppola etaH (120061) an d one from the sedentary 
survey by Giommi et al.l (|2005l ). 

The main selection criteria are the measured 
X-ray flux at 1 ke V and the distance of the ob- 
jects. According to IStecker et al.l (|1996l ). the syn- 
chrotron flux in the X-ray regime is connected to 
the flux in the VHE regime by 



VxFx ~ ^TeV-Frw > 



(1) 



assuming comparable synchrotron and Compton 
peak luminosities. Therefore objects with high 
X-ray fluxes are promising candidates for TeV 
emission. As the absorption of 7-rays within 
the extragalactic back ground light (EBL, see e.g. 
Kneiske fc Dole! 120101) is energy dependent, it is 
particularly important in the VHE regime to avoid 



strong attenuation of 7-r ays by limiting the red - 
shift range. According to Kneiske fc Dole! ( 2010h . 
at a redshift of z = 0.4, the expected cutoff en- 
ergy lies well above 200 GeV, allowing MAGIC to 
observe still with its highest sensitivity. There- 
fore all objects with a maximum redshift z = 0.4 
where considered. The energy threshold of the 
observations increases with the zenith distance 6. 
Accounting for this effect, the selection of sources 
with higher 9 (30° < 9 < 45°)during culmination 
should be limited to z < 0.15. The increasing ef- 
fect of EBL absorption should, however, imprint 
itself by a net steepening on the spectrum of the 
stacked excess. 

All criteria are described in detail below. They 
have been chosen to enhance the probability to 
detect the sources, hence we selected objects with 
high fluxes and inverse Compton peaks as well as 
allowing for the lowest possible energies to be mea- 
sured with MAGIC. 



Compared to lAlbert et al 



(|2008aD . the selec- 
tion criteria have been extended. The reason is 
the enhancement of the sample by taking a wider 
redshift or zenith distance range into account and 
including sources whose fitted synchrotron peak 
flux is high enough even if they show a lower X- 
ray flux level at 1 keV. The sample is divided into 
four parts: 

I. X-ray se l ected HB Ls obtained from 



Don ato et al.l (| 2 lh and lCostamante fc Ghisellini 
(|2002j): (i) redshift z < 0.4, (ii) X-ray flux 
-F^lkeV) > 2 /iJy, and (hi) zenith dis- 
tance 9 < 30° during culmination. As- 
suming the same luminosities at 1 keV as 
at 200 GeV (fol l owing the argumentation 
of IStecker et ail Il996l )). the X-ray flux 
i^(lkeV) = 2^iJy corresponds to a 7-ray 
flux at 200 GeV of - 4.8- 10~ 12 erg cm- 2 s" 1 . 
This criterion applies to 15 sources including 
nine sources already observed during cycle 1 
of regular MAGIC observations. The sources 
are listed in Table [T] 

• II. Two HBLs obtained from the same com- 
pilations taking a wider range in declina- 
tion and a lower maximum of the redshift 
into account: 1ES 0033+59.5 and RXS 
J1136. 5+6737. Selection criteria: (i) red- 
shift z < 0.15, (ii) X-ray flux F x {lkeV) > 
2 /iJy, and (hi) 9 < 45° during culmination. 
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• III. Intermed i ate B L Lacs taken from 
Nieppola et al. ( 20061 ) with high peak lu- 
minosities at the synchrotron peak. Se- 
lection criteria: (i) redshift z < 0.4, (ii) 
synchrotron peak frequency f p0 ak > 2 • 
10 15 Hz, (iii) flux at the peak F„ peak > 
10~ n crgcm~ 2 s _1 , and (iv) zenith angle 
6 < 30° during culmination. This is valid for 
three sources: B2 1215+30, PKS 1424+240, 
and B3 2247+381. All of them can also 



be found in iDonato et al.l ([20011 ) but with 
an X-ray flux at 1 keV below 2 fiJy. B2 
1215+30 is listed there as LBL. In return 
it is included in the TeV candidate list in 



Costamante fc Ghisellinil ((2002) 



• IV. One HBL from the sedentary survey 
(|Giommi et al.l 120051) with the same selec- 
tion criteria as applied for point one of the 
sample: 1RXS J044127.8+150455. 

As several other blazars fulfilling these selection 
criteria were already detected with MAGIC or 
other VHE instruments, a post-priori selection was 
done using only the objects which were not yet 
detected in the VHE regime in advance of the 
MAGIC observations leaving 21 objects as dis- 
cussed herein. All blazars in the MAGIC AGN 
observation program that fulfill these selection cri- 
teria either have been detected (or were known in 
advance) or they are listed here as non-detections. 

Tabled] lists all sources in the sample with rele- 
vant parameters. In case of multiple flux or spec- 
tral slope measurements the mean value is dis- 
played. 

3. OBSERVATIONS AND DATA ANAL- 
YSIS TECHNIQUE 

The observations presented here were carried 
out between December 2004 and April 2009 with 
a total amount of observation time of 490.0 hours. 
After quality selection (removing low quality data 
runs from the analysis) 394.1 hours were used for 
the analysis or 18.8 hours per source on average. 
The main reason to discard data from the analysis 
is a low event rate after image cleaning which is 
primarily influenced by the weather conditions. 

Most of the data were taken in wobble mode. 
In this mode the pointing position of the telescope 
is displaced by 0.4° from the source position. In 



order to get a well-balanced coverage inside the 
camera, the wobble position is changed regularly 
to the opposite (with respect to the source po- 
sition). Signal and background events are then 
determined from the same shower images with re- 
spect to the source position and to three symmet- 
ric OFF positions, respectively, all at the same 
distance to the camera center. Part of the data 
of RX J0319.8+1845, 2E 1415.6+2557 and RX 
J1725. 0+1152 were taken in ON mode where the 
pointing position of the telescope is centered on 
the object in the sky. For these observations ded- 
icated OFF observations have been used for the 
background estimation. 

The data were processed with the software 



package MARS (|Bretzl |2005h using an auto- 



mated analysis p i peline . Details can be f o und i n 
Bretz fc Waeneil 120031) iBretz fc Dornerl (|2008h . 



and I Albert et al.l (|2008b ). Furthermore, the ar 



rival time informat ion of neighborin g pixels was 



taken into account ( Aliu et al. 2009f) 



For the separation of signal and background 
events, dynamic cuts on the distribution of im- 
age parameters are applied. The image parame- 
ters are moments up to third orde r in the ligh t 
distribution of the shower images (Hillas 1985). 
The background suppressi on is done by mea ns of 
a parabolic cut in AREA (jRiegel et al.l 120051 ) and 
a cut in d 2 . The latter parameter is the squared 
angular distance between the source position and 
the reconstructed sho wer origin determined with a 
refined DISP method (jLessard et al.ll200ll ) taking 
into account the timing information of the show- 
ers. The z? 2 -cut used in this analysis is i? 2 < 0.0196 
which is a somewhat smaller value than usually 
used for the Crab Nebula, but provides a better 
background rejection for weak point sources. The 
chosen value for $ 2 corresponds to a signal region 
in the camera plane with a diameter of 2.8 cam- 
era pixels. The optical point spread function of 
the MAGIC telescope during the campaign was 
smaller than 16.0 mm corresponding to a diameter 
of 1.1 pixels, well within this area. A large sample 
of objects spanning a long time of observations has 
to be treated with a robust analysis. The usage 
of dynamic cuts provides such an analysis on the 
expense of sensitivity (cf. Section l4~Tj) . 

The statistical significance for any excess is 
calculated from the i9 2 distribution of signal 
and background events making use of Eq. 17 in 
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Table 1 
List of targets 



Object 


Season 


z 


logK) a 


F b 


F x c 


a x c 


Cat. d 


Sel. 












[MJy] 






crit. e 


lES 0033+595 


08/2006 - 


07/2008 


0.086' 


18.9 


2.0 


5.66 




C', N 


II 


1ES 0120+340 


08 - 


09/2005 


0.272 


18.3 


2.5 


4.34 


1.93 


C, D T , G, N 


I 


lES 0229+200 g 


08 - 


11/2006 


0.140 


19.5 


1.6 


2.88 




C*, N 


1,11 


RX J0319.8+1845 8 


12/2004 - 


01/2006 


0.190 


17.0 


0.4 


1.76 


2.07 


d\ g, n 


I 


lES 0323+022 


09 - 


12/2005 


0.147 


19.9 


6.3 


3.24 


2.46 


/") "T-v "I" X" "1 TV T 

C, D T , G, N 


1,11 


lES 0414+009° 


12/2005 - 


01/2006 


0.287 


20.7 


10.0 


5.00 


2.49 


C, D T , G, N 


I 


lRXS J044127.8+150455 


10 - 


12/2007 


0.109 






4.74 


2.10 


Gt 


IV 


lES 0647+250 


02 - 


03/2008 


0.203 f 


18.3 


3.2 


6.01 


2.47 


C T , D, N 


I 


lES 0806+524 g 


10 - 


12/2005 


0.138 


16.6 


1.6 


4.91 


2.93 


C, B\ N 


IJI 


lES 0927+500 


12/2005 - 


02/2006 


0.188 


21.1 


5.0 


4.00 


1.88 


Dt, G, N 


I 


lES I0ll+496 g 


03 - 


04/2006 


0.212 


16.7 


1.3 


2.15 


2.49 


C, Dt, N 


I 


lES 1028+5H 


03/2007- 


02/2008 


0.361 


18.6 


1.3 


4.42 


2.50 


C, Dt, G, N 


I 


RGB J1H7+202 


01/2007- 


03/2008 


0.140 






6.93 


1.90 


Ct, D, G 


1,11 


RX JH36.5+6737 




02/2007 


0.135 


17.6 


1.3 


3.17 


2.39 


C, Dt, G, N 


II 


B2 1215+30 


03/2007- 


03/2008 


0.237 


15.6 


1.3 


1.59 


2.65 


C, D, Nt h 


III 


2E 1415.6+2557 


04/2005 - 


04/2008 


0.237 


19.2 


3.2 


3.26 


2.25 


C, Dt, G, N 


I 


PKS 1424+240 S 


05/2006 - 


02/2007 


0.160 f 


15.7 


1.0 


1.37 


2.98 


D, Nt h 


III 


RX J1725.0+H52 


04/2005 - 


04/2009 


0.018 f 


15.8 


2.0 


3.60 


2.65 


C, Dt, N 


1,11 


lES 1727+502 


05/2006 - 


05/2007 


0.055 


17.4 


1.3 


3.36 


2.61 


C, Dt, N 


1,11 


1ES 1741+196 


07/2006 - 


04/2007 


0.083 


17.9 


1.0 


1.92' 


2.04 


C, Dt, N 


1,11 


B3 2247+381 


08 - 


09/2006 


0.119 


15.6 


1.0 


0.60 


2.51 


D, Nt h 


III 



a Fitted peak frequency from iNieppola et al.l (|2006| ) in units of log(Hz). 

b Flux at peak frequency extracted from INieppola et al. (2006) in units of 10 erg cm s . 
c Flux and photon spectral index at 1 keV. 

d Compilation wher e th e object appea r s (C: ICostamante fc Ghisellini (120021) : D: iDonato et all (l200lh : N: 
Nieppola et al. ( 20061 ); G: Giommi et al. ( 20051 )). The catalog from which the object was selected is marked 



with a dagger. 

°Selection criteria which are met by the object. 
'Tentative redshift. 

s Known VHE blazar (as of April 2010) due to a detection after the MAGIC observation period. 

h The objects chosen from INieppola et al.l ( 2006 ) are also listed in Donato et al. ( 2001 ). but with a X-ray 
flux lower than 2 /iJy. 



'Mean X-ray flux of multiple measurement in IDonato et al.l (|2001l ) below 2 /x Jy. 



Note. — List of objects in the sample of X-ray selected blazars with their observation time windows, redshifts 
and X-ray measurements. 
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Li fc Mai (119831) . 

Concerning the stacking method as described 
in Sections 14.41 and 1431 the -d 2 distributions have 
been summed up to retrieve the stacked d 2 signal 
plot. The differential energy spectrum is then cal- 
culated from all excess events using average values 
for the effective collection area and a Monte Carlo 
correction factor (spill-over correction), weighted 
each with the exposure time t, 



quality selection is t c 



19.2 h. The individual 



cxp ■ 



The same 

method has been applied to a data set of the Crab 
Nebula (cf. Section |4~T|) demonstrating its feasibil- 
ity. 

4. RESULTS OF THE MAGIC OBSER- 
VATIONS 

During the observation campaign no signifi- 
cant detection of any individual object could be 
achieved. The results can be found in Table [3] 
None of the objects showed flaring activity in 
the VHE band on a significant level on diurnal 
timescales within the observation time windows. 
Flaring activity is defined here as an offset of 3 
standard deviations from the mean measured 7- 
rate for each object. However, flux variations by a 
factor three would still prevent an individual ob- 
ject of the sample of being detected with high sig- 
nificance. In this Section we present the upper 
limits obtained for all 21 objects. 

Three of the objects were partially observed 
during an optical high state within a target of op- 
portunity campaign. The trigger criterion was an 
increase in the optical flux of the core of more 
than 50%. The objects are 1ES 0033+595, RGB 
J1117+202, and B2 1215+30. Significant activity 
or variability in the VHE 7-ray regime could not 
be detected. 

4.1. Crab Nebula Observations 

For a comparative analysis a sample of the Crab 
Nebula data has been used spanning a time range 
from Oct 2005 to Jan 2008. Three data sets 
have been chosen to account for the three differ- 
ent hardware conditions during the blazar observa- 
tions: 300 MHz readout system without and with 
optical splitters and 2 GHz readout system, later 
on referred to as 300 MHz, 300MHz O s and 2 GHz 
systems, respectively. The 9 distribution of the 
subsamples have been matched to the one of the 
blazar sample; the overall observation time after 



values as well as the combined result can be found 
in table [2] The energy spectrum can be fitted with 



a log p arabola (according to Eq. 2 in I Albert et al 
2008a ) accounting for the flattening of the spec- 



trum towards the inverse Compton peak: 



— = ( E 
dE ~ Io ' V300GcV 



\ [a+fclogio(E/300GeV)] 

J (2) 



with f = (5.37 + 0.11) ■ HT 10 TeV -1 cm" 2 s"\ 
a = -2.20 ± 0.05 and b = -0.11 ± 0.03. The z? 2 
distribution and the energy spectrum have been 
calculated in the very same way as for the blazar 
sample by stacking the three individual Crab Neb- 
ula samples. The integral flux above 150 GeV is 
determined to 

F>i50GeV = (2.81 + 0.05) • 10- 10 cm- 2 s- 1 . It will 
be used for comparison with the integral upper 
limits derived from the blazars. Figure [TJ displays 
the energy spectrum of the stacked excess of the 
Crab Nebula in comparison to the published spec- 
trum. The integral flux abo ve 150 GeV amounts t o 

"~ d2008ch . 



91 % of that determined in Albert et al 



A comparison to previous measurements of ex- 
periments like HEGRA, H.E.S.S or Whipple is 
difficult because of the higher energy threshold 
of these measurements (above 400 GeV). Due to 
the hardening of the Crab spectrum towards the 
peak below 100 GeV a simple extrapolation of the 
power-law spectra found there, overestimates the 
flux at 150 GeV leading to integ ral flux ratios of 
~ 70 - -80 % above 150 GeV (cf. lAharonian et al " 
2000;, l200d: lGrubell2007l) . 



The standard MAGIC integral sensitivity is 
- 1 .6 % of the Crab Nebula flux above 280 GeV for 
detecting a signal with 5 a in 50 hours (|Aliu et al. I 
2009). Including lower energies in the integral sen- 
sitivity determination, the value increases. The 
analysis presented in this work has an integral sen- 
sitivity above 150 GeV of 3.8 % of the Crab Nebula 
flux. This is mainly due to the long-term charac- 
teristics of the observations, because the analysis 
is aimed at a robust and conservative treatment of 
the data; in addition, also data before the instal- 
lation of the 2 GHz system are considered, where 
the standard MAGIC sensitivity above 280 GeV is 
also less with ^ 1.9 % of the Crab Nebula flux. 
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Table 2 

Observations of the Crab Nebula 



Season 


FADC 


^exp 




e 


Excess 


Backgr. 


Sign. 


-Ethr 




system 






[°] 


events 


events 


a 


[GeV] 


10/2005 - 03/2006 


300 MHz 


3.8 


6 


- 37 


967 


209 


36.0 


165 


09/2006 - 01/2007 


300MHz os 


8.1 


7 


- 43 


2086 


523 


51.0 


165 


02/2007 - 01/2008 


2 GHz 


7.3 


8 


- 30 


2133 


455 


53.5 


165 


Combined 




19.1 


6 


- 43 


5188 


1188 


82.2 


165 



Note. — Observations of the Crab Nebula used for a performance test of the stacking 
method and comparison to the flux upper limits of the blazars. The final spectrum (cf. 
Eq. [3]) is obtained as combination of all the subsamples. 
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Fig. 1. — Observations of three data sets of the Crab 
Nebula between October 2005 and January 2008. The red, 
blue and magenta colored data represent the data sets of 
the 300 MHz, 300MHz O s and the 2 GHz systems, respec- 
tively. The black curve shows the combined energy spec- 
trum obtained with the stac king method. For com parison 
the published spectrum from [Albert ct al. (2008^) is plot- 
ted as dashed green line. Note that the vertical error bars 
are hidden by the marks. 



4.2. Upper Limits 

The upper limits (U.L.) on the excess rates are 
calculated on a confi dence level of 3 a (9 9.7 %) us- 
ing the method from lRolke et al. ( 2005 ). Integral 
flux upper limits above a given energy are then 
calculated from them. The integral flux for each 
source is given above the energy threshold of the 
analysis, which is defined as the maximum of the 
differential distribution dN j dE vs E of simulated 
7-showers surviving all cuts. The integral fluxes 



are also compared to the integral flux of the Crab 
Nebula above the individual thresholds. 

The energy estimation for each source was done 
based on Monte Carlo simulated 7 events follow- 
ing a power-law distribution with T = —3.0 for 
a power law dN/dE oc E r . This was done in 
order to fit better the average spectral slope for 
the blazars in the VHE regime. For the integral 
upper limit calculation the same input spectrum 
(r = —3.0) was used. The resulting upper lim- 
its vary between 1.6% and 13.6% of the Crab 
Nebula flux above the individual energy thresh- 
old. The energy thresholds lie between 120 GeV 
and 230 GeV, due to differences in the 8 distribu- 
tions of the individual data samples. The Monte 
Carlo simulations have been chosen to match ex- 
actly the 8 distribution of each data sample. The 
results of the spectral analysis can be found in Ta- 
ble [31 too. 

Discovery of VHE 7-rays from RX J0319. 8+1845 
and 1ES 0806+524 has recently been reported 
by the VERITAS collaboration dAcciari et al 



2009; lOng k Fortinl l2009h . as well as from PKS 



1424+240 which was confirmed by the MAGIC 
collaboration in a campaign independent of the 



observations p resented here (jAcciari et al.1 l201dt 
Teshimal 120091 ). The measured VHE flux for 



the latter source was significantly higher than 
in previous observations with the MAGIC tele- 
scope. 1ES 0229+200 and 1ES 0414+009 have 
been detected by the H.E.S.S. tel escope ar- 



ray in 2006 (jAharonian et al.l l2007bT ) and 2009 
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Table 3 
Results of the analysis 



Object 


^exp 

[hi 

[ J 




9 

L J 


Excess 
events 


Backgr. 
events 


Scale 


Sign. 

(7 


rp a 
-C/thr 

[GeVl 


U.L. 
c.u. b 


U.L. 
f.u. c 


1ES 0033+595 


5.2 


31 


- 41 


60.0 


331.0 


0.33 


2.8 


170 


9.7 


2.4 


1ES 0120+340 


10.7 


6 


- 18 


20.7 


437.3 


0.33 


0.9 


120 


8.2 


3.1 


1ES 0229+200 


8.0 


8 


- 37 


55.0 


572.0 


0.33 


2.0 


120 


13.6 


5.1 


RX J0319. 8+1845 


11.2 


10 


- 31 


-23.4 


631.4 


0.59 


-0.7 


120 


1.6 


0.6 


1ES 0323+022 


11.4 


26 


- 46 


-45.3 


751.3 


0.33 


-1.5 


170 


6.9 


1.7 


1ES 0414+009 


18.2 


28 


- 36 


71.3 


1020.7 


0.33 


1.9 


170 


7.7 


1.9 


1RXS J044127.8+150455 


26.9 


13 


- 36 


18.3 


1825.7 


0.33 


0.4 


120 


3.2 


1.2 


1ES 0647+250 


29.2 


3 


- 32 


64.3 


1797.7 


0.33 


1.3 


120 


4.3 


1.6 


1ES 0806+524 


17.5 


24 


- 36 


17.0 


752.0 


0.33 


0.5 


140 


7.2 


2.2 


1ES 0927+500 


16.7 


21 


- 26 


28.3 


702.7 


0.33 


0.9 


140 


5.6 


1.7 


1ES 1011+496 


14.5 


21 


- 29 


89.0 


590.0 


0.33 


3.1 


140 


6.9 


2.1 


1ES 1028+511 


37.1 


22 


- 36 


65.7 


2312.3 


0.33 


1.2 


140 


3.3 


1.0 


RGB J1117+202 


14.9 


8 


- 38 


25.7 


804.3 


0.33 


0.8 


120 


5.3 


2.0 


RX J1136. 5+6737 


14.8 


39 


- 46 


22.7 


954.3 


0.33 


0.6 


230 


5.7 


0.9 


B2 1215+30 


16.1 


1 


- 41 


119.0 


995.0 


0.33 


3.2 


120 


9.3 


3.5 


2E 1415.6+2557 


57.4 


3 


- 36 


7.6 


3805.4 


0.54 


0.1 


120 


3.5 


1.3 


PKS 1424+240 


20.0 


5 


- 36 


51.7 


1210.3 


0.33 


1.3 


120 


8.2 


3.1 


RX J1725.0+1152 


32.0 


17 


- 35 


70.0 


1859.0 


0.38 


1.4 


140 


4.2 


1.3 


1ES 1727+502 


6.1 


21 


- 36 


31.0 


302.0 


0.33 


1.5 


140 


11.8 


3.6 


1ES 1741+196 


11.8 


9 


- 40 


98.7 


731.3 


0.33 


3.1 


120 


9.6 


3.6 


B3 2247+381 


8.3 


10 


- 36 


21.7 


490.3 


0.33 


0.8 


140 


5.2 


1.6 



a Pcak response energy for a power law spectrum E r with L = —3.0 

b Integral flux above E t hr given in units of the flux of the Crab Nebula (crab units, c.u.) 

c Integral flux above E t hr given in flux units f.u. = 10~ n cm~ 2 s _1 

Note. — Results of the analysis. The upper limits span a range of 1.6 - 13.6% of the Crab Nebula 
flux above the corresponding energy threshold. 



8 



( Hofmann fc Fegan 2009h . respectively. However, 
since the observations presented here were each 
performed in advance of the detections mentioned 
above, the inclusion of these sources in the stack- 
ing method is justified. The later detections show 
that the X-ray selection of possible targets is a 
reasonable approach. 

In order to compare the measured integral 
fluxes with the upper limits presented here they 
are extrapolated to the individual energy thresh- 
olds as reported in Table [3] In all cases except for 
PKS 1424+240 the upper limits are compatible 
with the extrapolated reported integral fluxes. 

4.3. Significance Distribution 

Taking a look at the calculated significances of 
the blazar sample it is evident that most of the in- 
dividual objects show positive values. Plotting the 
distribution of the significances, the mean value is 
not located at as expected for sky regions where 
no 7-rays are expected to originate. 

In Figure [2] the significance distribution for the 
blazar sample is shown together with the result of 
a cross-check as described below. As the num- 
ber of individual samples is different for both 
distributions they have been normalized to one. 
The blazar sample distribution has a mean value 
of 1.23±1.17 while the cross-check sample has 
— 0.08±0.85. This result can be expected due to 
the fact that our sample is biased by the selection 
toward potential VHE 7-ray emitters. 

In order to test if the positive signal in the 
blazar sample originates from a systematic effect 
of the observations or analysis chain, we also cross- 
checked this result with data sets obtained as OFF 
pointings associated to different ON source obser- 
vations not treated in this paper. These data sets 
were taken under similar conditions as the blazars 
covering the whole range of 9 of the blazar sample 
and processed with the very same analysis chain. 
The OFF observations were analyzed in wobble 
mode with respect to two fake source positions 
in the camera displaced by 0.4°from the camera 
center. Table 0] gives a list of these observations 
and results. Although the fit parameters of Gaus- 
sian fits to both distributions do not permit any 
conclusive statement, a Kolmogorov-Smirnov test 
of the compatibility of the blazar with the cross- 
check sample gives a probability of 1.56 %. For the 



Gaussian distributions the test returns 3.42 % and 
77.03 % for the compatibility of the blazar and the 
cross-check sample with the standard Gaussian, 
respectively. The cross-check sample is ~ 7 times 
smaller than the blazar sample, thus systematic ef- 
fects in the analysis can only largely be ruled out 
as possible explanation for the shift in the blazar 
distribution. 

4.4. Stacking Analysis 

Even if none of the sources was detected in 
a single observation, a cumulative signal search 
seems promising. For this reason the $ 2 -plots of 
the individual analyses have been stacked produc- 
ing one plot for the whole set containing 394.1 
hours of data (cf. Section [3]). Figure [3] shows the 
result, a significance of 4.9 standard deviations 
with 870 excess and 22876 background events. 
About 30% of the stacked excess comes from 
blazars now known as VHE 7-ray emitters. With- 
out these sources the stacked excess amounts to 
608 excess events with a significance of 3.8 a in- 
dicating that there are other emitters contained 
in the sample. Figure 2] underlines this finding. 
As expected, the stacked $ 2 -plot of the cross- 
check analysis containing no 7-signal gives a sig- 
nificance of —0.1 with —6 excess and 3009 back- 
ground events, the result is shown in Figure [3] as 
well. 

4.5. Energy Spectrum 

From the combined excess a differential energy 
spectrum can be calculated. The differential en- 
ergy spectrum dF/dE for one source is calculated 
binwise by dividing the product of the number of 
excess events N ea;Cj i and the spill over factor a^ by 
the product of effective collection area A e ff,i and 
exposure time t exp . In order to derive an energy 
spectrum of the stacked excess, the mean values 
of &t and Aoffj weighted with the observation time 
have to be taken: 



(at) 



En a «>™ ' ^' 



exp,n 



y^yf ^efF,i,n ' ^cxp,n 



zJn ^exp,n 



(3) 
(4) 



with n being the number of objects to be stacked 
and the energy bin i. The differential quotient 
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-t0.35 




Significance 



Fig. 2. — Significance distributions of the blazar (red, hatched up left to low right) and the cross-check sample (blue, hatched 
low left to up right). The different distributions are normalized to one, so the vertical axis gives the percentage of the whole 
blazar or cross-check sample, respectively. The blazar sample distribution has a mean value of 1.23±1.17 and the cross-check 
sample — 0.08±0.85. For comparison a Gaussian with mean value and standard deviation 1 is plotted as black curve. 



dNi/dE for each bin can then be calculated as 

dNj ^ En j^excM ■ (di) 

dE ~ E„ *exp,n • (^cff,i> • AEi U 

with the energy bin width AEj. The mean en- 
ergy spectrum in the observer's frame for all 21 
blazars considered in the stacking analysis can be 
well described by a power law, 

dN , , u 1 /£\^ m51 

(6) 

with E = 200GcV. The differential flux at 
200 GeV corresponds to 1.9% of the one of the 
Crab Nebula. The integral flux above 150 GeV 
is determined as F = 4.3 • 10~ 12 cm~ 2 s -1 corre- 
sponding to 1.5 % of the integral Crab Nebula flux 
above 150 GeV. 

On average, each blazar contributes with 
(2.1±0.3)/h excess events to the cumulative ex- 
cess as is illustrated in Figure 2] The objects are 
ordered in right ascension. 



In Figure [5] the measured spectrum is shown. 

5. DISCUSSION 

The positive mean significance distribution 
indicates that the X-ray selected blazars stud- 
ied here constitute a fairly representative sam- 
ple of generic VHE em i tters, as suggested by 
Costamante fc Ghisellinil (2002). The recent dis- 
coveries of individual blazars from the sample in- 
deed corroborate this finding. The next generation 
of Cherenkov experiments - MAGIC-II, H.E.S.S. 
2, and later on the Chere nkov Telescope Array 
(CTA, IWagner et all 120091) - will therefore have 
good chances to detect an increasing fraction of 
all known X-ray blazars. 

5.1. Gamma- ray Background 

At 200 GeV, the attenuation caused by the 

EBL is negligib le, according to the model of 

Kneiske fc Dole d2010h . so the calculation of 
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Table 4 
Data samples for cross-check 



Sample 




Season 


^cxp 


f\ 
U 




Excess 


Backgr. 


Sign. 








[h] 


[° 




events 


events 


a 


1 


06 


- 07/2006 


5.4 


34 - 


43 


-1.3 


335.3 


-0.1 


2 




07/2006 


3.1 


6 - 


29 


4.3 


107.7 


0.4 


3 




11/2006 


1.9 


37 - 


47 


19.0 


255.0 


1.0 


4 




01/2007 


3.3 


49 - 


56 


-24.7 


149.7 


-1.8 


5 




04/2007 


2.8 


11 - 


27 


-9.7 


139.7 


-0.7 


6 




05/2007 


1.3 


28 - 


37 


2.0 


76.0 


0.2 


7 




05/2007 


7.3 


29 - 


36 


-20.7 


356.7 


-1.0 


8 


01 


- 08/2008 


17.9 


22 - 


38 


7.0 


1041.0 


0.2 


9 


02 


- 04/2008 


9.3 


22 - 


26 


18.0 


548.0 


0.7 



Note. — Data samples used for the cross-check analysis. They were 
chosen to give a good coverage of the 9 distributions and the different 
night sky background conditions of the blazar sample. 



the broad-band spectral index ax_ 7 between 
1 keV and 200 GeV can be done with the ob- 
served VHE energy spectrum. The mean en- 
ergy flux at 200 GeV is calculated from the fit 
to 1.60 • 10 -12 ergcm~ 2 s _1 . This value is com- 
pared to the mean X-ray energy flux at 1 keV for 
all sources, weighted with their individual obser- 
vation time, which is 3.74 /xJy corresponding to a 
flux of 9.05-10- 12 ergcm- 2 s" 1 . The ratio of X-ray 
(IkeV) to 7-ray (200 GeV) flux is 



uF v (200 GeV) 



5.66 



(7) 



resulting in a broad-band spectral index ax 1 = 

1.09. 

The result suggests that during quiescence the X- 
ray luminosity is higher than the VHE 7-ray lu- 
minosity above 200 GeV. Here, we tacitly assume 
that the X-ray data, which are not contempo- 
raneous with the 7-ray data, are representative 
of baseline emission as well. Note, that the X- 
ray as well as the VHE data are an average over 
the whole blazar sample considered here and that 
flux variations commonly observed with the X-ray 
band do not influence ax-^ across eight orders 
of magnitude. A simple estimation of Aax_ 7 by 
inferring the error of the average value at 1 keV 



of the sample and the error of the energy spec- 
trum at 200 GeV results in Aa x -y = 0.04. With 
the newly found X-ray to 7-ray spectral index of 
ax-y = 1-09 one can infer the luminosity function 
of VHE blazars from their X-ray luminosity func- 
tion, avoiding the bias toward flares. Assuming 
equal X-ray and VHE 7-ray luminosities, HBLs 
already fail to e xplain the extragalactic diffu se 7- 
ray background (jKneiske fc Mannheim! 12008). 



5.2. Spectral energy distribution 

As no flaring activity has been seen on diurnal 
scales nor on longer time scales, the cumulative 
signal of the high-peaked blazars in this sample 
can be accounted as an upper limit on their base- 
line emission in VHE 7-rays, although variability 
on flux scales below the sensitivity limit of MAGIC 
may not be excluded. 

The SED for the blazar sample is determined by 
taking archival data in the radio and X-ray bands 
(1.4 GHz, 5 GHz, and IkeV) if available as well 
as contemporaneous optical data in the R-Band 
(640 nm) taken with the KVA telescope. The col- 
lected data are shown in Figure [6j In the VHE 
regime also the deabsorbed spectrum as already 
shown in Figure[5]is displayed. From the mean val- 
ues in optical and X-rays and the mean X-ray spec- 
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Blazar sample 
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Fig. 3. — ^-distribution of excess events for the stacked 
blazar sample (top) and the cross-check sample (bottom). 
The blazar sample shows a clear extension at low values 
with a significance of 4.9 standard deviations. 



• 600 




Fig. 4. — Excess events of the individual blazars vs. the 
overall exposure time. On average each blazar contributes 
with 2.1±0.3 excess events per hour. 
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tral index one can infer an average synchrotron 
peak energy of the sample below 1 keV. 

For a simple comparison the measured spec- 
tral energy distribution of the HBL 1ES 1959+650 
is drawn. 1ES 1959+650 is a well known VHE 
blazar which was observed in a historic low emis- 
sion state in a multiwa velength campaign in 2006 
(|Tagliaferri et al.1 120081 ) . The differential energy 



spectrum measured by MAGIC in the VHE regime 
follows a power law with a photon index T = 
-2.58 ± 0.18. The SED of 1ES 1959+650 can 
be well fitted with a one-zone synchrotron self- 
Compton model, which is also plotted in Figured] 
To guide the eye, the SED is also scaled down to 
the lowest energy bin of the VHE spectrum of the 
blazar sample. 



Fig. 5. — Differential energy spectrum obtained from the 
stacked source analysis. It is well described by a power law 
with index —3.16 ± 0.51. The integral flux above 150 GeV 
corresponds to 1.5% of the flux of the Crab Nebula. The 
spectrum of the Crab Nebula is shown as dashed gray line. 



6. CONCLUSIONS 

In the course of the MAGIC observational pro- 
gram during 2004 - 2009, a major part was spent 
on X-ray bright BL Lacertae objects. For 21 non- 
dctcctions upper limits on the integral flux ranging 
between 1.6 % and 13.6 % of the Crab Nebula flux 
could be determined. Applying a stacking method 
to the individual non-detections we found an av- 
erage VHE emission of the sample of X-ray se- 
lected blazars at the 4.9 a significance level above 
100 GeV. It turns out out that the mean VHE 7- 
ray flux is significantly lower than in archival X- 
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Fig. 6. — Spectral energy distribution vF v vs v for the blazar sample. Plotted are the measured fluxes in the radio (1.5 GHz 
and 5 GHz), optical (640 nm), and X-ray (IkeV) bands, together with their mean values. In the VHE regi me the observed 
spectr um of the stacked excess is shown. The black curve represents the model fit for 1ES 1959+650 taken from Taeliafcrri ct al. 
j200ST ) for comparison. The green curve is the same curve scaled down to match the first VHE flux point of the energy spectrum 
of the stacked excess. 



ray measurements. The two-point spectral index 
between IkeV and 200 GeV is 1.09±0.04. 
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